##########################################################
# John Henderson and Alex Theodoridis
# Replication Data for: "Seeing Spots", 
#  Forthcoming in Political Behavior, August 20, 2017
# 
##########################################################
#
#  figure_6.R 
#  -- file produces the estimates by discretion plot Figure 6 
#
##########################################################

rm(list=ls())
source('~/Dropbox/Seeing_Spots/replication/pre_data.R')     

ymat=cbind(video_skipped, replay, share, getlink, all_y,time_watched, total_time)   
 
newsint = as.numeric(newsint>2)

# pos v. neg & pov v. neg x interest          

set.seed(1005)   
pos_low_neg_high=pos_high_neg_low=list()
ovr=hvl=pvn=list() 
ivo=list()

for(k in 1:5){
	ovr[[k]]=pvn[[k]]=hvl[[k]]=matrix(NA,2000,2)
	pos_low_neg_high[[k]]=pos_high_neg_low[[k]]=matrix(NA,2000,2) 
	ivo[[k]]=matrix(NA,2000,2)
}

ym=as.matrix(cbind(all_y,video_skipped,replay,share,getlink))
#colnames(ymat)=c('all','skip','replay','share','getlink')

for(j in 1:2000){
	inds=sample(1:length(tr),size=length(tr),replace=T)
	omega=tr[inds]   
	ymat=ym[inds,]
	ominterest=newsint[inds] 
	ompid=pid_lean[inds]
	
	# pos_high_neg_low       	
	ind1 = which(str_sub(omega,1,3) == 'Pos' & ominterest == 1)
	ind2 = which(str_sub(omega,1,3) == 'Neg' & ominterest == 0)

	# pos_low_neg_high
	ind3 = which(str_sub(omega,1,3) == 'Pos' & ominterest == 0)
	ind4 = which(str_sub(omega,1,3) == 'Neg' & ominterest == 1)    
	
	# obama v. romney v. in v. out
	ind5 = which(ompid ==  1 & str_sub(omega,10,12) == "Oba")
	ind6 = which(ompid == -1 & str_sub(omega,10,12) == "Oba")
	ind7 = which(ompid == -1 & str_sub(omega,10,12) == "Rom")
	ind8 = which(ompid ==  1 & str_sub(omega,10,12) == "Rom")
	
	s1=sort(c(ind1,ind3)) # positive
	s2=sort(c(ind2,ind4)) # negative
	
	s3=sort(c(ind1,ind4)) # high
	s4=sort(c(ind2,ind3)) # low    
	
	s5=sort(c(ind5,ind6)) # obama
	s6=sort(c(ind7,ind8)) # romney  
	
	s7=sort(c(ind6,ind8)) # in
	s8=sort(c(ind5,ind7)) # out
	
	for(k in 1:5){
		pvn[[k]][j,1:2] = c(mean(ymat[s1,k],na.rm=T),mean(ymat[s2,k],na.rm=T))
		hvl[[k]][j,1:2] = c(mean(ymat[s3,k],na.rm=T),mean(ymat[s4,k],na.rm=T))
		  
		ovr[[k]][j,1:2] = c(mean(ymat[s5,k],na.rm=T),mean(ymat[s6,k],na.rm=T))		
		ivo[[k]][j,1:2] = c(mean(ymat[s7,k],na.rm=T),mean(ymat[s8,k],na.rm=T))  

		pos_high_neg_low[[k]][j,1:2] = c(mean(ymat[ind1,k],na.rm=T),mean(ymat[ind2,k],na.rm=T))
		pos_low_neg_high[[k]][j,1:2] = c(mean(ymat[ind3,k],na.rm=T),mean(ymat[ind4,k],na.rm=T))

	}  
}
    
colnames(ymat)=c('all','skip','replay','share','getlink')    
pvn[[1]]=pvn[[1]]/4
hvl[[1]]=hvl[[1]]/4  

ovr[[1]]=ovr[[1]]/4
ivo[[1]]=ivo[[1]]/4

pos_high_neg_low[[1]]=pos_high_neg_low[[1]]/4 
pos_low_neg_high[[1]]=pos_low_neg_high[[1]]/4 
      
ivoMU=c(mean((1-ivo[[2]][,1])-(1-ivo[[2]][,2])),
	    mean((ivo[[3]][,1])-(ivo[[3]][,2])),   
	    mean((ivo[[4]][,1])-(ivo[[4]][,2])),   
	    mean((ivo[[5]][,1])-(ivo[[5]][,2])))
ivoCI1=c(quantile((1-ivo[[2]][,1])-(1-ivo[[2]][,2]),prob=.975),
	     quantile((ivo[[3]][,1])-(ivo[[3]][,2]),prob=.975),   
	     quantile((ivo[[4]][,1])-(ivo[[4]][,2]),prob=.975),   
	     quantile((ivo[[5]][,1])-(ivo[[5]][,2]),prob=.975))
ivoCI0=c(quantile((1-ivo[[2]][,1])-(1-ivo[[2]][,2]),prob=.025),
	     quantile((ivo[[3]][,1])-(ivo[[3]][,2]),prob=.025),   
	     quantile((ivo[[4]][,1])-(ivo[[4]][,2]),prob=.025),   
	     quantile((ivo[[5]][,1])-(ivo[[5]][,2]),prob=.025))


pvnMU=c(mean((1-pvn[[2]][,2])-(1-pvn[[2]][,1])),
	    mean((pvn[[3]][,2])-(pvn[[3]][,1])),   
	    mean((pvn[[4]][,2])-(pvn[[4]][,1])),   
	    mean((pvn[[5]][,2])-(pvn[[5]][,1])))
pvnCI1=c(quantile((1-pvn[[2]][,2])-(1-pvn[[2]][,1]),prob=.975),
	     quantile((pvn[[3]][,2])-(pvn[[3]][,1]),prob=.975),   
	     quantile((pvn[[4]][,2])-(pvn[[4]][,1]),prob=.975),   
	     quantile((pvn[[5]][,2])-(pvn[[5]][,1]),prob=.975))
pvnCI0=c(quantile((1-pvn[[2]][,2])-(1-pvn[[2]][,1]),prob=.025),
	     quantile((pvn[[3]][,2])-(pvn[[3]][,1]),prob=.025),   
	     quantile((pvn[[4]][,2])-(pvn[[4]][,1]),prob=.025),   
	     quantile((pvn[[5]][,2])-(pvn[[5]][,1]),prob=.025))
         

md=rbind(pvnMU,ivoMU)
colnames(md)=c('Less Active','','','More Active')

# total means
total_mu=colMeans(ymat,na.rm=T)[2:5]
total_mu[1]=1-total_mu[1]

draw.lines=T
{
pdf(file='~/Dropbox/Seeing_Spots/replication/figures/discretion_barplot_US_lines.pdf',width=7,height=6)  
barplot(axes=F,as.matrix(md),beside=T,col=c('grey21','grey21'),density=c(22,500),ylim=c(-.115,.3),border=T,ylab='Difference in Selective Receipt',xlab='')
legend(x=1,y=c(.3),col=par('midnightblue'),legend=c('By Tone','By Source'),fill=c('grey21','grey21'),density=c(22,500),cex=1.1)
#text(x=c(2.1,5.1,8.1,11.1),y=c(-.065,-.065,-.065,-.065),labels=c('Not skip','Replay','Share','Get link'))
text(x=c(2.1,5.1,8.1,11.1),y=c(-.075,-.075,-.075,-.075),labels=c('Not skip','Replay','Share','Get link'))
axis(2)
#lines(x=c(1,12),y=c(-.1,-.1),lwd=2)
arrows(x0=c(.8),x1=12.3,y0=c(-.1),y1=-.1) 
arrows(x0=c(12.3),x1=.8,y0=c(-.1),y1=-.1)  
             
if(draw.lines==T){    
lines(x=c(1.5,1.5),y=c(pvnCI0[1],pvnCI1[1]),lty=2)
lines(x=c(1.5-.1,1.5+.1),y=c(pvnCI0[1],pvnCI0[1]),lty=1)
lines(x=c(1.5-.1,1.5+.1),y=c(pvnCI1[1],pvnCI1[1]),lty=1)

lines(x=c(2.5,2.5),y=c(ivoCI0[1],ivoCI1[1]),lty=1,col='white')
lines(x=c(2.5,2.5),y=c(ivoCI0[1],ivoCI1[1]),lty=2) 
lines(x=c(2.5-.1,2.5+.1),y=c(ivoCI0[1],ivoCI0[1]),lty=1,col='white')
lines(x=c(2.5-.1,2.5+.1),y=c(ivoCI1[1],ivoCI1[1]),lty=1)

lines(x=c(4.5,4.5),y=c(pvnCI0[2],pvnCI1[2]),lty=2)  
lines(x=c(4.5-.1,4.5+.1),y=c(pvnCI0[2],pvnCI0[2]),lty=1) 
lines(x=c(4.5-.1,4.5+.1),y=c(pvnCI1[2],pvnCI1[2]),lty=1) 

lines(x=c(5.5,5.5),y=c(ivoCI0[2],ivoCI1[2]),lty=1,col='white')
lines(x=c(5.5,5.5),y=c(ivoCI0[2],ivoCI1[2]),lty=2)      
lines(x=c(5.5-.1,5.5+.1),y=c(ivoCI0[2],ivoCI0[2]),lty=1,col='white')
lines(x=c(5.5-.1,5.5+.1),y=c(ivoCI1[2],ivoCI1[2]),lty=1)

lines(x=c(7.5,7.5),y=c(pvnCI0[3],pvnCI1[3]),lty=2)  
lines(x=c(7.5-.1,7.5+.1),y=c(pvnCI0[3],pvnCI0[3]),lty=1)
lines(x=c(7.5-.1,7.5+.1),y=c(pvnCI1[3],pvnCI1[3]),lty=1) 

lines(x=c(8.5,8.5),y=c(ivoCI0[3],ivoCI1[3]),lty=1,col='white')
lines(x=c(8.5,8.5),y=c(ivoCI0[3],ivoCI1[3]),lty=2)     
lines(x=c(8.5-.1,8.5+.1),y=c(ivoCI0[3],ivoCI0[3]),lty=1,col='white')
lines(x=c(8.5-.1,8.5+.1),y=c(ivoCI1[3],ivoCI1[3]),lty=1)                                               
                                                    
lines(x=c(10.5,10.5),y=c(pvnCI0[4],pvnCI1[4]),lty=2)      
lines(x=c(10.5-.1,10.5+.1),y=c(pvnCI0[4],pvnCI0[4]),lty=1)   
lines(x=c(10.5-.1,10.5+.1),y=c(pvnCI1[4],pvnCI1[4]),lty=1)   

lines(x=c(11.5,11.5),y=c(ivoCI0[4],ivoCI1[4]),lty=1,col='white')
lines(x=c(11.5,11.5),y=c(ivoCI0[4],ivoCI1[4]),lty=2)  
lines(x=c(11.5-.1,11.5+.1),y=c(ivoCI0[4],ivoCI0[4]),lty=1,col='white')
lines(x=c(11.5-.1,11.5+.1),y=c(ivoCI1[4],ivoCI1[4]),lty=1)   
}
dev.off()
}
  
# END figure_6.R